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1. Introduction 

Interacting quantum systems are incredibly complex, allowing a range of states that span 
an enormous Hilbert space. Some of these quantum states, such as squeezed and entangled 
states, can be utilized for quantum information [1,2] and precision measurement [3] purposes, 
potentially providing practical improvements over classical techniques. For instance, many 
common measurements based on quantum mechanical principles are usually limited to a much 
lower precision than allowed by the fundamental Heisenberg limit. The so-called standard 
quantum limit can be obtained by using the quantum analogue of classical states, such as 
the coherent state of light [4], to probe quantities ranging from positions and velocities, to 
magnetic fields and atom numbers. One can improve the precision by using entangled or 
squeezed states as inputs to interferometric measurement procedures [3, 5]. 

Experimentally realizing such quantum states presents many technical challenges. 
Macroscopic superpositions (or Schrodinger cat states) can in principle be used to perform 
Heisenberg-limited measurements [5], but in practice are extremely difficult to engineer. 
There has been much success in creating squeezed states of light [6] and the spin of atomic 
samples [7] that can be used to make measurements more accurate than the standard quantum 
limit. Interactions between samples of photons and/or atoms can create these squeezed states, 
but excess technical and thermal noise place practical limitations on the amount of squeezing 
generated. In order for quantum information and precision measurement techniques based on 
squeezed states to be practical, these excess fluctuations must be reduced. 
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Several methods of creating relative number squeezing and entanglement in ultra-cold 
gases have been proposed, and some realized, such as four- wave mixing [8, 9, 10] and 
the Kerr effect [11], molecular dissociation [12, 13], Josephson junction analogies [14, 15] 
and coupling with squeezed light [16]. Recently, Esteve et al. reported number difference 
squeezing and entanglement between multiple ultra-cold atomic clouds [15]. The experiment 
involved loading repulsive rubidium-87 atoms into two or more potential wells at low 
temperatures. The repulsion between the atoms reduces atomic bunching, and at zero 
temperature causes the fluctuations in the population differences between the wells to be 
smaller than the 'classical' binomial distribution. This experiment was the first to perform 
accurate atom counting within each well, while simultaneously having access to phase 
correlations between the wells through expansion and interference measurements. Combining 
these procedures, they observed relative number squeezing that could be used in principle to 
perform measurements at a precision 3.8 dB better than the standard quantum limit. 

One of the cornerstone features of this experiment is a technique employed to reduce 
the level of thermal fluctuations of the quantum state. Producing multiple condensates by 
evaporative cooling into multiple wells produced results close to the standard quantum limit. 
Esteve et al. found improved results by first cooling the atoms in a single trap, before 
modifying the potential to split the condensate between several wells. They postulated 
that this improvement is due to an adiabatic following of the quasiparticles representing the 
(quantized) fluctuations. Briefly, the energy of the excitations causing fluctuations is higher in 
the single condensate than after the splitting; therefore, the number of quasiparticle excitations 
at a given temperature is lower prior to splitting. As the potential is modified, the system is 
driven out of thermal equilibrium while the number of quasiparticles remains approximately 
fixed, and the level of fluctuations becomes lower than what would result from a thermal state 
at the minimum realizable temperature. In reference [15], the authors test this hypothesis by 
comparing with a simple, two-mode model, with reasonable agreement between the results. 

Here we present a critique of the two-mode model and the hypothesis of adiabatic 
passage. Going beyond the two mode model reveals regimes where two modes are sufficient, 
or when the approximation fails. Higher-order spatial modes contribute to the physics when 
the wells are not well separated, if the atom number is not large, or at higher energy or 
temperature scales. We employ the perturbative Bogoliubov method [17] taking into account 
the full multi-mode structure of the system and find the two-mode description is quite accurate 
for well-separated clouds. We compare the assumption of adiabatic passage to a model of 
full rethermalization, where the system remains in thermal equilibrium as the clouds are 
separated while the entropy of the isolated quantum system is constant. As the entropy of the 
total system depends on contributions from each spatial mode, the full multi-mode statistics 
need to be calculated for this comparison. We find that only the adiabatic model produces 
the dramatic improvement in squeezing observed in the experiment. Finally, we perform 
truncated Wigner [18] simulations of the dynamics and find they agree well with the adiabatic 
model for a sufficiently large system. Thus we can conclude that thermalization between 
the quasiparticle modes in this system is slow, and adiabatic passage is a good model for 
experiments performed on these timescales. 

2. System model 

The experiment begins by evaporatively cooling a cloud of *'^Rb atoms in an optical double 
well potential. The double-well potential is created by superimposing an optical lattice onto 
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Figure 1. A schematic of the double well potential (dashed red) and ground-state mean- 
field density (solid black) when the chemical potential /i is close to the barrier height. The 
ground state can be interpreted as a superposition of the atoms being in the left and right wells, 
sketched by the respective dotted-line "wave-functions". 



a harmonic dipole trap [19, 15]. The potential in the region of the condensate is given by 

V{x, y, z) = ^mulix^ + y^) + ^mw^z^ + cos(2A:lz), (1) 

where the atomic mass of '^^Rb is m w 1.44 x 10~^^ kg. The trap is cylindrically symmetric 
and elongated in the direction of the lattice, uj^ < uj±. The trap has radial trapping frequency 
uj± w 27r X 425 Hz and, for the two-well experiment, a longitudinal trapping frequency of 
uiz 2tt X 60 Hz. The optical lattice creates a barrier between two (or more) low energy 
regions, or wells, such as that depicted in figure 1 and is generated by interfering two phase- 
locked lasers at an angle producing a lattice spacing of a « 5.7 fim, where the lattice wave- 
vector is = 7r/a. The peak-to-peak depth Vl begins at 27rfi, x 430 Hz before being ramped 
up. 

The atoms are described by the Hamiltonian, 

' -ft' 2 . ... Uo 



H = ip^ 



-V^ + V{r,t) + ^7/-^^ 



d-'r. (2) 



2m ' ' ' 2 

The interaction constant is Uq — ATrti^as/m, where the s-wave scattering length is as ~ 5.39 
nm. 

We are interested in the low temperature limit, at or near the ground state of the system. 
To the lowest order of approximation, the solution is given by the ground-state mean-field 
wave-function ipQ, which we find by integrating the imaginary-time Gross-Pitaevskii equation, 

_^d^ = ^VVo + (Vir) + C/olVoH^o - /^V-o, (3) 
OT Zm 

until the state reaches equilibrium. The chemical potential fj, is adjusted to give the required 
atom number A'o = / \ipo\'^ d'^r = 1600. If/iis close to the barrier height, the solutionis that 
of two weakly-coupled condensates, sketched in figure 1 . 

The peak-to-peak depth of the optical lattice begins at 27r/ix 430 Hz. As the gas is cooled, 
a Bose-Einstein condensate forms in the trap with the potential barrier between the wells 
significantly smaller than the chemical potential /i — thus at this stage it can be considered 
a single condensate. The evaporative cooling method is unable to reduce the temperature T 
significantly below the level of the chemical potential [15], so fc^T ~ fi. 

The potential barrier is slowly ramped up by increasing the intensity of the lasers. In 
the experiment, this was achieved adiabatically (i.e. slow enough that additional heating was 
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not observed) at a rate where the peak-to-peak lattice potential Vl increases by approximately 
27r?i X 2 Hz per ms. The BEC is coherently split in two, whereupon the atoms configure 
themselves in a low energy state possessing relative number squeezing. 

The experiment was then able to directly detect the atoms in situ, and count the atoms in 
each well, Ni to a greater accuracy than the shot-noise limit (i.e. AiV^ < ^/Nl). Defining the 
relative number squeezing, or normalized variance, as 

e = Var[iVi-7V2]/(iVi+7V2), (4) 

their results indicate squeezing of up to -6.6 dB, or a reduction in the number difference 
variance by a factor of 4.5 compared to the binomial distribution expected from ideal, 
uncorrected condensates. To quantatively determine the number variance we must go 
beyond the Gross-Pitaevskii description, which does not include either quantum or thermal 
fluctuations in the cloud. 



3. Bogoliubov analysis 

The Bogoliubov approach is a perturbative expansion valid in the weakly-interacting or 
large atom number limit [17, 20]. In either case, the ground state of the system is close 
to a coherent state given by the Gross-Pitaevskii equation. The Bogoliubov approach treats 

fluctuations about the mean-field perturbatively, by writing ^(r, t) = [?/'o(r)+(5?/'(r, t)]e~*^*, 

and assuming Sijj is 'small' compared to ?/'o. The linearized solutions for Si) arising from 
equation (2) obey 

JV7 + [/o^2^V, (5) 

One can diagonalize the above linear equation (which couples 5~t}) with ^V'S into their 
respective eigenstates. The annihilation operators for the eigenstates have the form 

= e*"'*/'"' j u*{r)5^{r) + v,{r)5il:\r)d'^r. (6) 

The eigenstates with eigenvalue ti are defined by the functions Ui{r) and Vi{r). Combining 
the above yields the Bogoliubov-de Gennes (BdG) equations 

" /:gp + C/o|V'oP -c/oV-o 

where for brevity we have defined the Gross-Pitaevskii operator £gp — ^^^^ + ^(^) + 
t^l^oP ~ We implement the procedure presented in reference [20] to reliably solve 
equation (7). 

To efficiently solve the BdG equations, it is important to minimize the size of our 
three-dimensional spatial basis. Here we implement the harmonic oscillator basis [21] 
as it effectively represents the trapped BEC in the lowest energy states while allowing 
efficient calculations of the matrix elements between the basis states (f>i{r), for example 
Vij = J (l)*{r)^ cos{2kLz)(f>j{r) d^r and Uij ~ J \-4'Q{r)\(l>*{r) (l>j{r) d^r can he found 
accurately using the Gaussian quadrature sum rules [22, 23]. We further optimize the 
calculations by taking advantage of x, y and z reflection symmetries, allowing us to reduce 
the computation cost by roughly a factor of 8^. The total system is truncated to the 2413 states 
having energy less than QOhuJz- 

We display the condensate density and lowest energy quasiparticle mode solutions in 
figure 2 for lattice depths of 430 Hz and 1650 Hz. This lowest energy mode can be 
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Figure 2. Results of the Bogoliubov-de Gennes approach at a lattice depth of (a,b,c) 430 Hz 
and (d,e,f) 1650 Hz. (a,d) Density of the mean field |i/'o(r)p through the x-z plane. As the 
lattice is ramped the condensate splits into two well-defined wells. (b,c,e,f) The solutions for 
the lowest energy excitation, (b,d) ui(r) and (c,e) vi(r), also through the x-z plane. This 
excitation removes particles from one well and transfers them to the other. 
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Figure 3. (a) Energies of the lowest quasiparticle excitations as a function of lattice potential 
depth Vl ■ The lowest energy state is indicated and bold, (b) The first-order contributions to 
the number difference variance Ai, as given by equations (9) and (10). The lowest energy state 
contributes most strongly to the number difference variance. 



interpreted as transferring particles from one well to the other and corresponds to Josephson- 
style oscillations in population. For deep lattices the lowest energy mode has a significant 
overlap with the antisymmetrized ground state wave-function; that is u{r) and v{r) are 
approximately proportional to tpoif) x P{z), where P{z) is —1 for z < and +1 otherwise 
(i.e. P{z) = 8(z) — Q{—z), where 8 is the Heaviside step function). In figure 3 we see that 
the energy of this lowest mode decreases as the lattice depth is increased, while the energy of 
the other modes increase. For Vl ^ ^nh x 1500 Hz, the energy of the second excited mode 
is several times that of the lowest excitations, and at low temperatures one would expect the 
majority of quasiparticles in the lowest mode. 

3.1. Calculating correlations 

Based on the solutions to the Bogoliubov-de Gennes equations, one can calculate the variance 
of the atom number difference between the two wells as a function of the total atom number, 
temperatures and the potential. We begin by defining the regions of interest around each well, 
r^i and ^2- For the double well we choose 51i to be the region with z > and il2 to be the 
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region with 2 < 0. The quantum operator describing the total number of atoms in the ith 
region is — ip'' {r)il} (r) (Pr. Expanding the field operator about the mean field using 

the solutions to the Bogoliubov-de Gennes equations, ■0(r) = ipoir) ~ H'^ii''')' 

and inserting into the above gives 

N., = \Mr)\' + J2Mr) {u*{r)b] - v*{r)b^) 
+ Y,^*{r)uk{r)b% - u*{r)vk{r)b]bl 

jk 

- v*{r)uk{r)b^b^ + v*{r)vk{r)bp^ d'r. (8) 

The system is symmetric under the transformation z ^ —z, and so we expect the wells 
to be evenly balanced, {N-^ — N2) = 0. We take an ansatz for the many-body density 
matrix by assuming the Bogoliubov modes are in thermal (chaotic) states, with population 
(SJ&j) = [exp(ei/fcBT) — 1]^^. The variance of the number difference is therefore 

j 

+ J2{\^^A'' + C^,D,,){b^b^){b;b]) 

+ E 1' + ^^^^^^ (^M) (9) 

where we used the symmetry of the system to simplify the expression, and defined the 
following integrals for brevity: 



A, = J Piz){roir)u,{r)-Mr)v*ir))d^r, (10) 

B^j = J P{z)u*{r)v*{r)d'r, (11) 

= j P{z)u*{r)u,{r)d^r, (12) 

D,, = [ P{z)v,{r)v*{r)d^r, (13) 



remembering P{z) is —1 for 2 < and +1 otherwise. We note that some the terms in equation 
(9) go beyond the accuracy of the first-order Bogoliubov method used to derive them, and 
many authors therefore choose to disregard these terms. Here these terms are retained, without 
loss of accuracy, to facilitate a better comparison with the truncated Wigner simulations in 
section 4, where the full density matrix resulting from thermally-occupied Bogoliubov modes 
is implemented as initial conditions. 

For large A'o and small fluctuations, the dominant contributions to the number variance 
come from the first sum in equation (9), and in particular the lowest energy (Josephson) mode 
has the largest contribution. In figure 3 (b), we see that for deep lattices |Ai p is much larger 
than those corresponding to other modes, and simultaneously the quasiparticle population in 
the lowest energy mode is expected to be greatest. Together these might justify the two mode 



Multimode analysis of non-classical correlations in double well BECs 



1 



Multimode, r = 
Multimode, r= 11 nK 
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Figure 4. The number difference squeezing ^ [equation (4)] as a function of barrier height for 
the full multi-mode analysis (black lines) and a simplified, two-mode model involving just the 
lowest Bogoliubov excitation (red lines). The dashed lines are the results at zero temperature, 
while the solid lines are for T = 11 nK. The models agree best for large lattice heights, or 
well-separated clouds. 
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Figure 5. A contour plot of the number difference squeezing g as a function of temperature 
and barrier height in the double well system. Low temperatures and large barriers produce the 
least relative number fluctuations. 



approximation for lattice intensities Vl ^ 27r/i x 1500 Hz, whereas higher energy excitations 
are important for weaker lattices. 

This would indeed be the case for systems of larger particle number. However, in this 
system the additional contributions to equation (9) that a two mode model would fail to 
capture are significant. Even for well separated wells {Vl ^ 27r/i x 1500 Hz) these extra 
terms contribute approximately one third of the total variance at zero temperature. In this 
regime, we expect a small loss of precision in our calculations due to the fact that those 
terms proportional to Bij, Cij and Dik go beyond the accuracy of the first-order Bogoliubov 
analysis used, and may be poorly approximated. 

In figure 4 we compare the results of the full multi-mode analysis presented here, and 
a simplified model including just the lowest Bogoliubov excitation. This simplified model is 
somewhat akin to a two-mode approximation, and we see that the two models show reasonable 
agreement for large barrier heights. In this regime, the two atomic clouds are well-separated 
and the two-mode approximation is expected to be accurate. 

In figure 5 we plot the number squeezing ^ [equation (4)] as a function of barrier height 
and temperature. Clearly, best squeezing is obtained at lower temperatures and with larger 
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barrier heights. 

3.2. Adiabatic Passage 

We now compare our results to the double-well experiment performed at Universitat 
Heidelberg and examine the hypothesis of adiabatic passage presented in reference [15]. 
The experiment found improved number squeezing when slowly raising the barrier after 
evaporation compared to performing evaporative cooling with the barrier at fixed height. If 
the lattice ramp is slow compared to the energy gaps in the Bogoliubov spectrum, ei, 62 — ei, 
et cetera, one might expect the population in each quasiparticle mode to be fixed during the 
evolution. If this were true, then the system would no longer be in thermal equilibrium, as 
the ratio of the excitation energies changes during the ramp (see figure 3). The ramp ends 
with fewer quasiparticle excitations in the lowest mode than one would expect at realizable 
temperatures. 

However, for sufficiently slow ramps one might expect the system to rethermalize, so 
that the population in each Bogoliubov changes to conform to a new global temperature as the 
energy levels move, and we also consider this possibility. The evolution of the potential, 
and therefore the quasiparticle energies, is slow compared to the gaps in the Bogoliubov 
spectrum, but it is unclear whether the evolution is adiabatic with respect to the full many-body 
Hamiltonian. The Bogoliubov description ignores interactions between the quasiparticles that 
can lead to a redistribution of excitations and a return to thermal equilibrium. 

An isolated system undergoing quasi-static evolution will have constant entropy. The 
entropy of mode i in a thermal state, in terms of the mean occupation, is [24] 



The entropy of the total system is S* = J^i ^i- the lattice height, and therefore excitation 
spectrum, varries, we calculate the new temperature required to keep the system in thermal 
equilibrium with fixed entropy. As many modes make significant contributions to the total 
entropy, and entropy is exchanged between modes, this analysis is not possible with a two- 
mode model. 

We compare the two models for thermalization as the lattice potential is increased 
in figure 6. In all cases we begin in thermal equilibrium with a lattice strength of 430 Hz at 
a temperature of ksT « /i/6, chosen to roughly correspond to the experimentally observed 
fluctuations. 

As can be seen from figure 6, the adiabatic process predicts a level of fluctuations 
significantly lower than the fully thermalizing model. In a realistic experiment one might 
expect the entropy to increase slightly, resulting in even greater fluctuations. As the 
experiment observes strongly improved number difference squeezing as the lattice is ramped, 
one can rule out the possibility of rethermalization on the timescales of the experiment. On 
the other hand, the assumption of adiabatic passage agrees well with the experimental results, 
and we conclude that the mechanism of effective cooling proposed in [15] is correct. 

4. Dynamical simulation 

In this section we perform a direct dynamical simulation using the truncated Wigner [18, 
25, 10, 21] method and compare the results to perfect adiabatic passage. As pointed 
out in the previous section, the number fluctuations of the double-well system realized 
by [15] are sensitive to second-order interactions between the Unearized BogoUubov modes. 




(14) 
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Figure 6. A compaiison of different thermalization models for the relative number squeezing 
as the lattice height is slowly ramped from 430 Hz to 1650 Hz (indicated by the vertical dotted 
lines). In all cases, at a lattice depth of 430 Hz we take a thermal distribution with T = fj,/6. 
The dashed red line is the squeezing where we assume the state remains in thermal equilibrium 
with constant entropy. The solid black line is for the adiabatic model where the number of 
quasiparticle excitations in each mode does not change. The adiabatic model agrees with the 
experimental observation of significantly decreased number fluctuations when slowly ramping 
up the lattice. 



For the experimental parameters, the truncated Wigner simulations show dynamics that 
is inconsistent with the linearized Bogohubov-de Gennes solutions - in particular, the 
quasiparticle populations and number difference variance are not constant in time even for 
a fixed Hamiltonian. The Bogoliubov fluctuations are too large to be considered perturbative, 
and the linearized approximation breaks down. 

Nevertheless, we can probe adiabatic passage dynamically by investigating a slightly 

different system where the perturbations S-ip are small compared to i/jq, and linearization is 
accurate. To minimize the changes to the system we investigate a system of more atoms with 
smaller interactions, such that the product NqUq is unchanged. In principle, such a change 
could be realized with a Feshbach resonance [17], although this would be challenging. 

The truncated Wigner method is an approximate phase space method that takes advantage 
of the Wigner representation of a quantum field. Truncating the higher-order derivatives from 
the equation of motion for the Wigner function results in a classical Liouville equation that 
can be efficiently sampled stochastically. Thus we reduced the problem of solving a non- 
linear equation for a quantum field i^{r) to that of the non-linear motion of an ensemble of 
classical fields ip{r). The resulting trajectories obey the (real-time) Gross-Pitaevskii equation, 
where the initial conditions are sampled from the initial Wigner distribution. Expectation 
values of symmetrically ordered operators are equal to the ensemble average of their classical 
counterparts; however care must be taken when using a finite, non-uniform basis [21]. 

The system we simulate has 1.6 x 10^ atoms and a scattering length Gs = 5.39 x 10^^^ 
m, and begins in thermal equilibrium at Vl = ^irh x 430 Hz with a temperature T « /i/6. 
This system has 100 times more atoms than the experiment, but identical mean-field and 
Bogoliubov-de Gennes solutions. The initial state is populated by a coherent state condensate 
surrounded by Bogoliubov modes in thermal (chaotic) states with populations given by the 
Bose-Einstein statistics [17]. The Wigner distribution for such a state is Gaussian and 
therefore straightforward to randomly sample. We evolve 4000 independent trajectories with 
the Gross-Pitaevskii equation, while the lattice is ramped up at a rate 27rfi x 2 Hz/ms. As we 
use the non-uniform Hermite-Gauss basis for our simulation, care must be taken to extract the 
correct quantum expectation values from the averaged data. For arbitrary basis {(piif)} such 
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Figure 7. A comparison of the truncated Wigner dynamics (points), a model of complete 
adiabatic following (black solid line), and the fully thermalized model (red dashed line) for 
the system with 1.6 X 10^ atoms and a small scattering length. The dynamics agree with the 
model of adiabatic passage. 



that 'ip{r) — "Y^ Ci4'i {f) ™d matrix P defined by 

= J cl,*{r)P{r)^,{r)d'r, (15) 
we find the number difference 

(N^-N^)^ J 4;^r)P{r)tP{r) = ctPc - Trace [P] /2, (16) 
and the number difference squared 

{{N^ - N^)^) = (ctPc- Trace[F]/2)^ - Trace[p2]/4. (17) 

These relations allow us to efficiently calculate the number difference variance without 
transforming to the spatial basis, while simplifying the symmetric Wigner corrections for 
a non-uniform basis [21]. 

Our dynamical results are plotted in figure 7, along with the predictions for perfect 
adiabatic following and full thermalization obtained by the same procedure as in the previous 
section. We see that the number difference variance agrees with adiabatic model, validating 
the hypothesis of adiabatic following on this timescale. 



5. Systems with different geometry 

It is worth discussing changes that occur in multi-well systems of different geometry. First, 
double well traps can be created in a variety of ways, such as the fully optical approach 
analysed here [15] or with magnetic traps on atom chips [26]. The latter setup results in two 
parallel, elongated condensates (with uj^ ^ and ojz). In such a system, one would expect 
long wavelength excitations along the transverse dimension x, some of which will have energy 
less than the Josephson mode. Similarly, one expects multiple Josephson-type modes within 
a small energy band, coupling the typical double-well excitation with low lying transverse 
modes (the limit of two infinite condensates can be found in [27]). In this case, a multimode 
description such as is presented here is essential to analysing these systems. It is unclear, 
however, if such a geometry is advantageous for creating sub-shot noise number correlations. 
The increased density of states about the Josephson mode may make adiabatic following 
more difficult. Furthermore, a clear definition of the relative phase of the condensates, as 
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investigated in reference [15], may be problematic as the phase may vary along the long axis 
of the system. 

In reference [15] the authors also investigated systems of multiple wells. They presented 
improved results for the inner pair of six coupled wells, which was theoretically investigated 
using a six-mode Bogoliubov theory. We performed a multi-mode analysis of this system, 
which involved an increased computational difficulty. In the multi-mode picture, we observed 
a band of five Josephson-type states of decreasing energy as the barrier height is increased 
— one mode for each inter-well coupling — while the remaining excitations had increasing 
energy. Thus, we can conclude that the adiabatic following technique should also be benificial 
in this arrangement, as was indeed verified experimentally [15]. We were unfortunately unable 
to repeat the full analysis for that system because of the increased numberical difficulty. 

6. Conclusion 

We have performed a multi-mode Bogoliubov analysis of coupled Bose-Einstein condensate 
systems to find the low temperature physics and calculate the population fluctuations between 
the two wells. At very low temperatures, repulsive interactions induce relative number 
squeezing and even entanglement between the wells. We present a critical analysis of the 
adiabatic passage hypothesis employed in reference [15] as the optical lattice separating 
the wells is ramped up. Our results indicate that adiabatic evolution of the quasiparticles 
is possible, causing a decrease in number fluctuations as the lattice height is increased. 
We conclude that the isolated system can be driven out of thermal equilibrium, to reduce 
fluctuations to below what is possible at thermal equihbrium. It remains to be seen if this 
procedure can be improved to reach the quantum regime, where the vacuum fluctuations 
dominate the physics. 

There are several Umitations to the Bogoliubov approximation; we are unable to 
investigate either the low atom number or high temperature limits. The transition between 
anti-bunching due to repulsion and boson bunching at higher temperatures could be 
investigated by non-perturbative methods. Quantum Monte Carlo techniques [28] provide 
the most accurate and reliable results for high temperature bosonic systems. The Positive-P 
phase space method can be implemented to find thermal statistics of a Bose gas [29]. The 
Popov approximation [17] is similar to the approach used here, but considers the interactions 
between the excited modes and back-action on the condensate, and may be an appropriate 
method for exploring this regime. Related to these, Sinatra et al. [30] have proposed a method 
for sampling the statistics of a Bose gas described by a Gaussian density matrix. 
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